subroutine free_energy(energy)
use type
use var
implicit none
integer :: i,j,k,label,step
real(dp) :: LdG = 0.d0, Elastic = 0.d0, surface0 = 0.d0, surface1 = 0.d0
real(dp), intent(out) :: energy(5) 

LdG = 0.d0
Elastic = 0.d0
surface0 = 0.d0
surface1 = 0.d0

  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           
           if (bulk(i,j,k) == .true.) then

              LdG = LdG + LdG_Energy(q(:,label))
              Elastic = Elastic + Elastic_Energy(label)
              
           elseif (boundary(i,j,k) == .true..and.infinite==.false.) then
              
              surface0 = surface0 + surface0_Energy(label)
           
           elseif (boundary_NP(i,j,k) == .true..and.NP_infinite==.false.) then
              
              surface1 = surface1 + surface1_Energy(label)
           endif

      enddo
    enddo
  enddo
  
  energy(1) =  LdG * vd
  energy(2) =  Elastic * vd
  if (nondegenerate) then
      energy(3) =  0.5d0*gama(1)*surface0*ad0
      energy(4) =  0.5d0*gama(2)*surface1*ad1
  elseif (degenerate) then
      energy(3) =  gama(1)*surface0*ad0
      energy(4) =  gama(2)*surface1*ad1
  endif
  energy(5) =  sum(energy(1:4))

return
contains

Function LdG_Energy(q_tem)
use type
use var
implicit none
real(dp), dimension(6) :: q_tem
real(dp) :: tr1, tr2, tr3
real(dp) :: aux1, a1, a4
real(dp) ::  LdG_Energy

tr1 = sum(fac*q_tem*q_tem)

tr2 = tr1 * tr1

aux1 = 3.d0/2.d0
a1 = -q_tem(1) - q_tem(4)
a4 = q_tem(1) - q_tem(4)

tr3 = 0.75d0*a1**3 - 3.d0*a1*q_tem(2)**2 + aux1*a1*q_tem(3)**2 + aux1*q_tem(3)**2*a4 &
       - 0.5d0*aux1*a1*a4**2 + 6.d0*q_tem(2)*q_tem(3)*q_tem(5) + aux1*a1*q_tem(5)**2 - aux1*a4*q_tem(5)**2
LdG_Energy = 0.5d0*( 1.d0 - ULdG/3.d0 )*tr1 - ULdG/3.d0*tr3 + ULdG/4.d0*tr2
return
end function LdG_Energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Function  Elastic_Energy(label_c)
use type
use var
use fix
implicit none
integer :: ii
integer :: x, y, left, right, label_l, label_r, label_c
real(dp), dimension(18) :: grad
real(dp) :: Elastic_Energy

Do x = 1 , 3
   ii = (x-1)*2 + 1
   label_l = neighbors(ii,label_c)
   label_r = neighbors(ii+1,label_c)
   grad(1+(x-1)*6:x*6) = (q(:,label_r)-q(:,label_l))*df2(x)
   grad(1+(x-1)*6:x*6) = grad(1+(x-1)*6:x*6)*grad(1+(x-1)*6:x*6)*fac
enddo

!grad = grad*grad*fac
Elastic_Energy = 0.5d0*sum(grad)

return
end function Elastic_Energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Function  surface0_Energy(label_c)
use type
use var
implicit none
integer ::i,indx, label_c
real(dp) :: rapini, surface0_Energy, galota
real(dp), dimension(6) :: Qin

indx = lb2index(label_c)
if (nondegenerate) then
   rapini = 0.d0
   rapini = sum (fac*(q(:,label_c)-qb0(:,indx))*(q(:,label_c)-qb0(:,indx)))
   surface0_Energy = rapini
elseif (degenerate) then
   galota = 0
   call surface(q(:,label_c),normal0(:,indx),Qin)
   galota = sum(fac*Qin*Qin)
  surface0_Energy = galota
endif
end function surface0_Energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Function  surface1_Energy(label_c)
use type
use var
implicit none
integer ::i,indx, label_c
real(dp) :: rapini, surface1_Energy, galota
real(dp), dimension(6) :: Qin

indx = lb2index(label_c)
if (NP_nondegenerate) then
   rapini = 0.d0
   rapini = sum (fac*(q(:,label_c)-qb1(:,indx))*(q(:,label_c)-qb1(:,indx)))
   surface1_Energy = rapini
elseif (NP_degenerate) then
   galota = 0
   call surface(q(:,label_c),normal1(:,indx),Qin)
   galota = sum(fac*Qin*Qin)
   surface1_Energy = galota
endif
end function surface1_Energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine surface(Qin,nu,Qout)
use var
use type
implicit none
real(dp), dimension(6) ,intent(in)::Qin
real(dp), dimension(3) ,intent(in)::nu
real(dp), dimension(6) ,intent(out)::Qout
real(dp), dimension(3,3) :: Qm, P, Qpp, Qtt
integer :: i,j,k,l
Qpp =0.d0;Qtt=0.d0
Qm(1,1)=Qin(1);Qm(1,2)=Qin(2);Qm(1,3)=Qin(3);Qm(2,1)=Qin(2);Qm(2,2)=Qin(4)
Qm(2,3)=Qin(5);Qm(3,1)=Qin(3);Qm(3,2)=Qin(5);Qm(3,3)=Qin(6)

Do i=1,3
   Do j=1,3
      if (i==j) Qtt(i,j) = Qm(i,j) + Sbulk/3.d0
      if (i/=j) Qtt(i,j) = Qm(i,j) 
      
      if (i==j) P(i,j) = 1 - nu(i)*nu(j)
      if (i/=j) P(i,j) = - nu(i)*nu(j)
   enddo
enddo

Do i=1,3
   Do j=1,3
      Do k=1,3
         Do l=1,3
            Qpp(i,j) = Qpp(i,j) + P(i,k) * Qtt(k,l) * P(l,j)
         enddo
      enddo
  enddo
enddo

Qout(1)=Qtt(1,1)-Qpp(1,1)
Qout(2)=Qtt(1,2)-Qpp(1,2)
Qout(3)=Qtt(1,3)-Qpp(1,3)
Qout(4)=Qtt(2,2)-Qpp(2,2)
Qout(5)=Qtt(2,3)-Qpp(2,3)
Qout(6)=Qtt(3,3)-Qpp(3,3)

endsubroutine surface



End subroutine  free_energy


